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By  F.  Ham  and  Y.-N.  Young 


1.  Motivation  and  objectives 

Simulations  of  flows  involving  free  surfaces  have  become  ubiquitous  in  the  literature. 
The  evolution  of  both  simulation  methods  and  computer  speed  have  allowed  the  inves- 
tigation of  many  problems  of  increasing  complexity  and  engineering  relevance.  For  2D 
simulations,  it  is  reasonably  common  to  use  grid  resolutions  of  512  x 512  or  larger.  See 
for  example  the  2D  planar  breaking  wave  simulations  of  Chen  et  al.  (1999),  or  the  ax- 
isymmetric  drop  breakup  simulations  of  Han  and  Tryggvason  (1999a,  1999b).  Many  such 
free  surface  problems  are  fundamentally  three  dimensional  and  the  absence  of  three- 
dimensional  effects  such  as  the  pinch-off  of  stretched  filaments  by  surface  tension  makes 
it  difficult  to  make  valid  physical  conclusions  or  use  2D  results  to  develop  models.  Achiev- 
ing similar  grid  resolution  in  3D,  however,  is  outside  the  computational  reach  of  most 
research  centers. 

For  many  problems  of  interest,  a fine  grid  is  required  in  only  a relatively  small  portion 
of  the  domain  (e.g.  around  regions  of  high  surface  curvature)  and  much  coarser  grids  are 
sufficient  to  capture  the  smoothly  evolving  velocity  fields  far  from  the  surface.  It  is  this 
type  of  problem  that  motivates  the  present  Cartesian  adaptive  method,  where  the  finest 
resolution  in  the  neighbourhood  of  the  free  surface  may  be  of  order  5123  (sometimes 
called  the  effective  resolution),  but  the  coarsening  away  from  the  free  surface  yields  an 
actual  number  of  control  volumes  that  is  less  than  using  the  fine  grid  everywhere  by  at 
least  an  order  of  magnitude. 

A number  of  adaptive  Cartesian  methods  for  freesurface  calculations  have  been  pro- 
posed in  the  literature.  For  example,  Haj-Hariri  et  al.  (1997)  used  Cartesian  local  grid 
refinement  to  investigate  thermocapillary  motion  of  3D  drops.  Their  method  used  an 
octree-based  structure  to  store  the  grid  heirarchy.  More  recently  Sussman  et  al.  (1999) 
proposed  an  adaptive  levelset  method  based  on  the  embedded  grid  approach,  where  the 
global  computational  domain  is  divided  into  uniformly  refined  rectangular  patches.  This 
approach  has  the  advantage  of  allowing  a traditional  staggered  or  semi-staggered  struc- 
tured discretization  throughout  most  of  the  domain,  with  special  interpolations  required 
only  at  the  interface  between  grid  levels.  The  embedded  grid  approach  may  not,  however, 
lead  to  the  smallest  Cartesian  grids  possible  because  of  the  restriction  that  refined  regions 
be  rectangular  and  isotropic.  Additionally,  embedded  grid  approaches  can  be  difficult  to 
parallelize  and  load  balance. 

In  the  present  contribution  we  develop  a level  set  method  based  on  local  anisotropic 
Cartesian  adaptation  as  described  in  Ham  et  al.  (2002).  Such  an  approach  should  allow 
for  the  smallest  possible  Cartesian  grid  capable  of  resolving  a given  flow.  The  remainder 
of  the  paper  is  organized  as  follows.  In  section  2 the  level  set  formulation  for  free  surface 
calculations  is  presented  and  its  strengths  and  weaknesses  relative  to  the  other  free  surface 
methods  reviewed.  In  section  3 the  collocated  numerical  method  is  described.  In  section 
4 the  method  is  validated  by  solving  the  2D  and  3D  drop  oscilation  problem.  In  section 
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5 we  present  some  results  from  more  complex  cases  including  the  3D  drop  breakup  in  an 
impulsively  accelerated  free  stream,  and  the  3D  immiscible  Rayleigh-Taylor  instability. 
Conclusions  are  given  in  section  6. 


2.  Mathematical  Formulation 

For  recent  reviews  of  the  level  set  method  and  its  many  applications,  the  interested 
reader  is  referred  to  Osher  and  Fedkiw  (2001)  or  Sethian  (2001).  Here  we  briefly  present 
its  formulation  for  incompressible  2-fluid  problems. 

The  incompressible,  immiscible,  two-fluid  system  is  treated  as  a single  fluid  with  strong 
variations  in  density  and  viscosity  in  the  neighborhood  of  the  interface.  The  continuity 
and  momentum  equations  for  such  a flow  can  be  written  in  conservative  form  as: 


dp 

dt 
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dpui 

dt 


dpuj 


= 0, 
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dxj 
dpujUi 
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dp  dra  ..  _ 
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(2.1) 
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where  Uj  is  the  fluid  velocity,  p the  fluid  density,  p the  pressure,  Ty  the  viscous  stress 
tensor,  gl  the  acceleration  due  to  gravity,  a the  surface  tension  coefficient,  k the  local  free 
surface  curvature,  S the  Dirac  delta  function  evaluated  based  on  d the  normal  distance 
to  the  surface,  and  the  unit  normal  at  the  free  surface. 

It  is  well  known  that  solving  equations  2.1  and  2.2  directly  in  the  presence  of  the 
discontinuities  at  the  interface  will  lead  to  the  excessive  smearing  of  the  interface  over 
long  time  integration,  and/or  numerical  “ringing”  in  the  region  of  the  discontinuities  due 
to  dispersive  errors  in  the  numerical  discretization.  In  the  present  work,  we  capture  the 
interface  as  the  zero-level  isosurface  of  a higher-order  function  p,  the  level  set  function 
(effectively  an  additional  transported  scalar).  On  either  side  of  the  interface,  the  sign  of 
the  level  set  function  can  be  used  to  determine  which  fluid  is  present.  This  allows  the 
fluid  properties  to  be  written  as  a function  of  the  level  set  only,  specifically  p = p{P)- 
This  allows  the  continuity  equation  to  be  written: 
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For  the  case  of  constant  density  except  in  the  neighborhood  of  the  zero  level  set, 
which  represents  the  fluid  interface,  the  solution  of  the  continuity  equation  can  thus  be 
decomposed  into  the  solution  of  the  following: 
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Eq.  2.4  is  the  standard  evolution  equation  for  the  level  set  function,  and  eq.  2.5  is  the 
incompressible  continuity  equation.  Note  that  the  level  set  equation  need  only  be  solved 
near  the  interface,  which  we  have  defined  as  the  zero  level  set.  Away  from  the  interface, 
it  is  common  to  make  p equal  to  a signed  distance  function. 

In  a similar  way,  the  application  of  chain  rule  to  the  momentum  equation  allows  the 
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Figure  1.  Spatial  location  of  variables  for  collocated  discretization. 


decomposition  of  the  momentum  equation  into  the  solution  of  the  level  set  equation,  eq. 
2.4,  and  the  solution  of  the  following  momentum  equation: 


dui  : dujUi 


dxj 


^Zl+^+p9,+,’KSmn, 


(2.6) 


Equations  2.4  - 2.6  comprise  our  governing  system  for  the  present  work.  Note  that 
the  level  set  formulation  results  in  a system  of  governing  equations  that  can  no  longer 
be  written  in  conservative  form.  Thus,  when  the  finite  volume  method  is  applied  to  this 
system,  we  cannot  expect  to  achieve  discrete  conservation  of  mass  and  momentum  (in 
the  region  of  the  interface).  It  is  the  hope  of  the  level  set  formulation  that  these  errors 
in  conservation  are  mitigated  by  the  more  accurate  capturing  of  the  interface  that  is 
possible  with  the  smoothly- varying  level  set  function. 


3.  Numerical  Method 

The  governing  equations  are  solved  on  a Cartesian  adaptive  grid.  The  grid  arrangement 
and  adaptation  are  described  more  completely  in  Ham  et  al.  (2002),  and  here  we  only 
discuss  details  specific  to  the  simulation  of  multiphase  flows  using  level  set  methods  on 
such  a grid. 

Figure  1 shows  the  spatial  arrangement  of  variables  on  the  Cartesian  adaptive  grid. 
All  variables  are  stored  at  the  control  volume  (CV)  centers  with  the  exception  of  a 
face-normal  velocity,  located  at  the  face  centers,  and  used  to  enforce  the  divergence- 
free  constraint  in  each  time  step.  The  variables  are  staggered  in  time  so  that  they  are 
located  most  conveniently  for  the  time  advancement  scheme.  Denoting  the  time  level  by 
a superscript  index,  the  velocities  are  located  at  time  level  tn  and  tn+1,  and  pressure, 
density,  viscosity,  and  the  level  set  at  time  levels  tn~i  and  tn+i . The  semi-descretization 
of  the  governing  equations  in  each  time  step  is  then  as  follows. 


Step  1.  Advance  the  level  set: 

cf>n+i 


<t>n~i 


+ 


“?5Z5-(*"_i  + *"+i) 


At 


0 


(3.1) 
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Here  the  advantage  of  time-staggering  is  evident  - i.e.  we  know  the  velocity  uf  from 
the  previous  time  step,  so  this  equation  is  decoupled  from  the  other  equations,  and  can 
be  advanced  on  its  own  to  get  the  level  set  value  at  the  mid-point  of  the  new  time  level, 
<f>n+i.  Here  we  have  used  an  implicit  Crank-Nicholson  time  advancement  scheme,  which 
requires  an  iterative  solution  method.  In  practice  the  hyperbolic  system  is  not  stiff,  and 
can  be  quickly  converged  by  a simple  iterative  scheme  such  as  Gauss-Seidel.  The  spatial 
derivatives  required  in  eq.  3.1  are  approximated  by  a 5t,l-order  WENO  scheme  (Jiang  & 
Peng  2000).  While  effective  at  recovering  the  viscosity  solution  to  the  level  set  equation, 
the  WENO  discretization  requires  that  the  grid  be  locally  refined  in  a uniform  way  in  a 
band  about  the  zero  level  set,  and  thus  does  not  take  full  advantage  of  the  flexibility  in 
the  adaptive  grid. 

Step  2.  Reinitialize  the  level  set  by  solving  to  steady  state: 


j£  = sgn  (</>)(  1-IWI).  (3.2) 

The  reinitialization  step  is  performed  every  time  step,  and  ensures  that  the  level  set  is  a 
signed  distance  function  away  from  <j>  = 0,  without  (theoretically)  changing  the  location  of 
the  zero-level  set.  Spatial  derivatives  required  in  eq.  3.2  are  once  again  approximated  using 
a 5t,l-order  WENO  scheme  (Jiang  & Peng  2000).  For  this  equation  we  use  explicit  third 
order  TVD  Runge-Kutta  method  for  time  integration.  In  practice,  it  is  not  neccessary 
to  solve  to  steady  state,  and  we  have  found  that  5 iterations  with  a pseudo-time  step  of 
At  = Ax/ 4 is  sufficient  when  performed  every  computational  time  step. 

Step  3.  Calculate  the  new  density  and  viscosity  at  the  midpoint  of  the  time  step: 


Pn+i  = PlH  (<T+  = ) +P2(l-H  (r+i)) 
//"+*  - mH  (<t>n+*)  +M2  (l  (<£n+5)) 


(3.3) 

(3.4) 


where  pi,p\  are  the  constant  properties  of  the  fluid  occupying  the  region  <f>  > 0,  and 
p2i  P2  are  the  constant  properties  of  the  fluid  occupying  the  region  <f>  < 0,  and  H is  a 
smoothed  Heaviside  step  function  (Sussman  et  al.  1994). 


Step  4.  Project  the  momentum  equations: 

The  remaining  steps  are  a variant  of  the  collocated  fractional  step  method  described, 
for  example,  by  Kim  and  Choi  (2000).  First,  calculate  a projected  velocity  field  Ui  using 
the  momentum  equations: 


*■<  1_  ( dp_n  * , pn+i 

At  pn+h  < 


(3.5) 


where  Ri  contains  all  other  terms  in  the  momentum  equation,  eq  2.6.  We  discretize  both 
the  convective  and  viscous  terms  implicitly  using  second-order  symmetric  discretizations, 
and  the  surface  tension  explicitly  based  on  the  known  level  set  at  the  midpoint  of  the 
current  time  step,  <j)n+% . Some  care  must  be  taken  in  the  discretization  of  the  surface 
tension  terms  when  treated  as  source  terms  in  the  momentum  equations  of  a collocated 
scheme,  as  described  below. 
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Step  5.  Subtract  the  old  pressure  gradient: 

u*n+l=~n+l  AtlSP_n 
1 p”+i  Sxi 

Step  6.  Interpolate  the  starred  velocity  field  to  the  faces: 


U}n+1  =u*n+lf  - At 


R^1  RnfH 

pn+i  n+}. 


where  ()  is  a second-order  interpolation  operator  that  yields  a face-normal  component 
from  two  CV-centered  vectors. 

Step  7.  Solve  the  variable  coefficient  Poisson  equation  for  the  new  pressure: 

5Er'Af=E  -h  r*A'-  <3-8> 

/ / Pf 

Step  8.  Update  the  face  velocities  to  the  new  divergence-free  field: 

ur'-Uf  i ,39, 

A t ~ »+§  dn  [ 1 

pf 

Step  9.  Reconstruct  the  pressure  gradient  at  the  CV  centers: 

_L_dP_n+i  _ _±_dpn+*Xi  (310) 

pn+i  dXi  "+*  dn  ^ 1 

H Pf 

where  ()  ’ represents  a reconstruction  operator.  In  this  case  we  use  a face-area  weighted 
least  squares  reconstruction. 

Step  10.  Update  the  CV  velocities: 


<+1-<  i_dp_n+*  ( . 

At  pn+k  dXi  } 

By  adding  eqs.  3.5,  3.6,  and  3.11,  it  is  clear  that  a second-order  time  advancement  of 
the  momentum  equation  is  recovered. 

A critical  difference  between  the  present  formulation  and  the  formulation  of  Kim  and 
Choi  is  in  the  calculation  of  the  starred  face-normal  velocities  (eq.  3.7).  Kim  and  Choi 
assume  that: 


Rn+l/2f  J^+V2 


p?1/2 ' 
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This  is  an  0( Ax2)  approximation,  seemingly  consistent  with  the  overall  accuracy  of  the 
method,  and  significantly  simplifies  the  calculation  of  the  Poisson  equation  source  term. 
In  the  present  investigation,  however,  it  was  found  that  when  surface  tension  forces  were 
introduced  in  the  region  of  the  zero  level  set,  this  approximation  could  lead  to  large  non- 
physical oscillations  in  the  CV-centered  velocity  field.  To  solve  this  problem,  the  surface 
tension  forces  must  be  calculated  at  the  faces,  and  then  averaged  to  the  CV  centers,  i.e.: 


P PS 


(3.13) 


With  this  calculation  of  the  surface  tension  forces,  we  can  no  longer  make  the  assump- 
tion of  eq.  3.12,  and  the  additional  terms  must  be  included  in  the  calculation  of  t/j?  and 
thus  in  the  Poisson  equation  source  term. 


4.  Validation 

As  validation  cases  for  the  method,  we  solve  the  2D  column  and  3D  drop  oscillation 
problems. 


4.1.  2D  column  oscillation 

Prom  Lamb  (1932),  for  a 2D  column  of  liquid  in  the  limit  of  zero  viscosity  perturbed 
according  to: 


r = a (1  4-  ecos  ( nO ))  (4.1) 

where  a is  the  mean  radius,  e a small  perturbation,  and  n the  mode,  the  oscillation 
frequency  will  be: 


u>2  = n (n2  — l)  -^r.  (4.2) 

paA 

Figure  2(a)  shows  the  initial  condition  for  a drop  oscillation  calculation  for  mode 
n = 5.  Here  we  have  exaggerated  the  perturbation  amplitude  to  e = 0.15  for  illustrative 
purposes.  In  the  actual  computations,  e = 0.02  was  used  for  the  initial  condition.  Figure 
2(b)  compares  the  calculated  results  to  the  theory  (period  T = r/w),  with  generally 

good  agreement.  For  these  calculations,  we  used  periodic  boundary  conditions  in  a square 
domain  of  size  L = 4a,  and  set  the  density  of  the  surrounding  fluid  to  p2  — O.OOlp  to 
minimize  its  effect  on  the  oscillations.  The  fluid  viscosity  was  adjusted  to  keep  the  decay 
in  amplitude  of  the  oscillations  to  less  that  5%  per  cycle. 

4.2.  3D  drop  oscillation 
For  a 3D  drop  perturbed  according  to 


r = a (1  + eSn ) 

where  the  surface  harmonics  of  order  n,  Sn,  are  defined: 


(4.3) 


(4.4) 
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mode  n 


(a) 


(b) 


FIGURE  2.  a)  Sample  grid  and  initial  zero- level  set  for  2D  drop  oscillation  problem  with  n = 5. 
b)  Comparison  of  computed  and  theoretical  2D  drop  oscillation  periods  for  modes  n = 2 through 
5. 


52=iVf(W0_1) 

(4.5) 

S3  = ^\j^~  (5 cos3 9 — 3 cosO') 

(4.6) 

S4  = — (35 cosi6  — 30 cosd  + 3) 

(4.7) 

S3  — (63  cos46  — 70  cos26  + 15) 

(4.8) 

• is  given  by: 

u>2  = n (n  — 1)  (n  + 2) 

paA 

(4.9) 

Figure  3 (a)  shows  the  drop  surface  for  the  mode  n = 5.  Figure  3 (b)  compares  the 
calculated  results  to  the  theory,  also  with  good  agreement. 


5.  Results 

In  the  following  subsections  we  present  some  results  from  simulations  of  more  com- 
plex flows.  For  the  purposes  of  this  paper  describing  the  numerical  method,  these  cases 
are  meant  to  simply  illustrate  the  potential  of  the  method.  Consequently  we  do  not 
present  any  analysis  of  the  associated  flow  physics,  which  will  be  the  subject  of  future 
investigations. 
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(a)  (b) 

Figure  3.  a)  Initial  zero-level  set  for  3D  drop  oscillation  problem  with  n = 5.  b)  Comparison 
of  computed  and  theoretical  3D  drop  oscillation  periods  for  modes  n = 2 through  5. 

5.1.  Secondary  breakup  of  spherical  drops 

The  secondary  breakup  of  drops  is  an  important  problem  in  spray  atomization,  and  the 
subject  of  the  2D  axisymmetric  numerical  investigation  of  Han  and  Tryggvason  (1999a, 
1999b).  We  used  the  present  Cartesian  adaptive  method  to  calculate  3D  drop  breakup  in 
an  impulsively  accelerated  free  stream.  Figure  4 shows  the  calculated  free  surface  for  a 
case  run  with  equivalent  resolution  of  128  x 128  x 128  with  a domain  size  of  5 D x 5 D x 5 D. 

5.2.  Rayleigh- Taylor  Problem 

The  Rayleigh-Taylor  instability  refers  to  the  instability  of  the  plane  interface  between 
two  fluids  of  different  density  superimposed  one  over  the  other  and  subject  to  gravity. 
When  the  upper  fluid  has  a density  greater  than  the  lower  fluid,  the  interface  can  be 
unstable  to  small  perturbations  whose  amplification  is  well  described  by  linear  theory 
with  dependence  on  the  density  ratio,  gravity,  surface  tension  coefficient,  and  viscos- 
ity (Chandrasekhar  1961).  As  the  amplification  continues  and  the  problem  enters  the 
nonlinear  regime,  the  fluid  mixing  can  become  chaotic,  characterized  by  bubbles  of  lighter 
fluid  rising  into  the  heavier  fluid,  and  spikes  of  heavier  fluid  falling  into  the  lighter  fluid, 
with  regions  of  high  vorticity  near  the  spike/bubble  interface.  The  resulting  mixing  zone 
is  known  experimentally  (Schneider  et  al.  1998)  to  broaden  in  a way  that  depends  lin- 
early on  the  gravitational  acceleration  g and  the  Atwood  number  A = (pi  — p2)/{p\  +P2), 
and  quadratically  on  the  time,  i.e. 


hbts  — Ub.sAgt  , (5.1) 

where  hb  or  h3  represents  the  penetration  length  for  bubbles  or  spikes  respectively,  and 
a has  been  introduced  as  a constant  of  proportionality,  sometimes  called  the  acceleration 
constant. 

Recently,  the  miscible  Rayleigh-Taylor  problem  has  been  investigated  numerically  by 
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(a)  0.5  (b)  1.5  (c)  2.5  (d)  3.5  (e)  4.5 


(f)  5.5  (g)  6.5  (h)  8.5 

Figure  4.  Breakup  of  an  initially  spherical  drop  in  an  impulsively  accelerated  field.  The 
number  below  each  figure  indicates  the  non-dimensional  time  tU/D.  Free  stream  veloc- 
ity is  from  top-left-back  to  lower-right-ffont.  Properties  for  this  simulation  are  as  follows: 
Re  = PgDU/ng  = 1000,  We  = pgU2D/a  = 20,  Oh  = pi/Vpi^D  = 0.007. 

Young  et  al.  (2001),  who  calculate  at  sa  0.03.  Experimental  results  for  immiscible  fluids 
yield  higher  values  of  a\,  « 0.05  — 0.07.  A recent  review  of  experimental  and  simulation 
results  for  the  immiscible  case  is  given  by  Glimm  et  al.  (2001). 

Figure  5 presents  some  results  from  the  application  of  the  present  Cartesian  adaptive 
method  to  this  problem.  These  computations  were  performed  in  a 2 x 2 x 2 box  with 
the  interface  perturbed  randomly  with  a maximum  amplitude  of  0.01.  Properties  of  the 
simulation  were  selected  to  give  a wavelength  of  maximum  instability  of  A*  = 0.35.  These 


Figure  5.  Evolution  of  the  interface  for  the  3D  Rayleigh-Taylor  problem  with  multi-mode 
perturbation:  top  panels  - plan  view  looking  down  on  interface;  bottom  panels  - elevation  view 
from  side. 

preliminary  simulations  were  made  with  the  relatively  coarse  equivalent  resolution  in  the 
region  of  the  zero-level  set  of  64  x 64  x 64,  or  about  10  CVs  per  wavelength  for  the  most 
unstable  mode.  Figure  5 illustrates  the  evolution  of  the  interface  at  3 different  times  for 
the  Atwood  number  A = 0.98. 

6.  Conclusions 

A simulation  tool  that  integrates  Cartesian  adaptive  grids  with  the  level  set  method 
has  been  developed.  The  method  as  described  is  particularly  suited  to  problems  where 
the  smallest  scales  are  associated  with  the  free  surface  motions,  and  significant  savings 
can  thus  be  realized  by  coarsening  the  grid  away  from  the  surface.  The  method  has  been 
validated  and  its  potential  demonstrated  by  solving  several  test  problems,  including  2D 
and  3D  drop  oscillations,  drop  breakup  in  an  impulsively  accelerated  free  stream,  and 
the  immiscible  3D  Rayleigh-Taylor  problem. 

Our  ongoing  work  is  focused  on  changing  the  level  set  solver  to  the  particle  level  set 
method  of  Enright  et  al.  (2002) . The  improved  conservation  properties  of  this  method 
may  allow  us  to  move  away  from  the  higher-order  WENO  discretizations  presently  used 
in  the  level  set  advancement  and  reinitialization  equations.  With  lower-order  methods, 
adaptation  can  be  added  in  the  region  of  the  zero-level  set,  resulting  in  a further  signifi- 
cant reduction  in  computational  cost. 
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